##############################################################################

library(data.table)
library(argparser)
library(dplyr)

##############################################################################

argp <- arg_parser("Qc Somatic variant")
argp <- add_argument(argp, "--Tumor" , help="")
argp <- add_argument(argp, "--Normal" , help="")
argp <- add_argument(argp, "--out_file" , help="")
argp <- add_argument(argp, "--skip_line" , help="")
argp <- add_argument(argp, "--maf_file" , help="")
argp <- add_argument(argp, "--qc_file" , help="")
argp <- add_argument(argp, "--stand_list_file" , help="")
argp <- add_argument(argp, "--out_record_file" , help="")

argv <- parse_args(argp)

Tumor <- argv$Tumor
Normal <- argv$Normal
skip_line <- as.numeric(argv$skip_line)
maf_file <- argv$maf_file
qc_file <- argv$qc_file
stand_list_file <- argv$stand_list_file
out_record_file <- argv$out_record_file
out_file <- argv$out_file

if(1!=1){

  Tumor <- "WGC079284"
  Normal <- "WGC079302"
  skip_line <- 3439
  maf_file <- paste0("/public/home/xxf2019/20220205_lungSomatic/results/maf/",Tumor , "_" , Normal,".filterStep1.maf")
  qc_file <- paste0("/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/",Tumor , "_" , Normal,"_MutQc.annotationRegion.csv")
  stand_list_file <- "/public/home/xxf2019/20220205_lungSomatic/scripts/mutect2_Scripts/QC_Standard.list"
  out_record_file <- paste0("/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/",Tumor , "_" , Normal,"_MutQc.RecordQcNum.csv")
  out_file <- paste0("/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/",Tumor , "_" , Normal,"_MutQc.KeepMut.csv")

}

##############################################################################

dat_maf <- fread(maf_file  , skip = skip_line , quote = "")
dat_qc <- fread(qc_file)
dat_stand <- fread(stand_list_file)

##############################################################################
## 统一INDEL的格式
dat_qc$Chromosome <- sapply( strsplit(dat_qc$POS , ":") , "[" , 1)
dat_qc$Start_Position <- sapply( strsplit(dat_qc$POS , ":") , "[" , 2)
dat_qc$Tumor_Seq_Allele1 <- sapply( strsplit(dat_qc$POS , ":") , "[" , 3)
dat_qc$Tumor_Seq_Allele2 <- sapply( strsplit(dat_qc$POS , ":") , "[" , 4)
dat_qc <- subset(  dat_qc , POS!="" )
dat_qc$mut_length <- abs(nchar(dat_qc$Tumor_Seq_Allele1) - nchar(dat_qc$Tumor_Seq_Allele2))

## INDEL
result <- c()
for( i in 1:nrow(dat_qc)){
  tmp <- dat_qc[i,]
  chr <- tmp$Chromosome
  pos <- as.numeric(tmp$Start_Position)
  ref_length <- nchar(tmp$Tumor_Seq_Allele1)
  alt_length <- nchar(tmp$Tumor_Seq_Allele2)

  if( ref_length < alt_length ){
    ## INS
    ref <- "-"
    alt <- substr( tmp$Tumor_Seq_Allele2 , 2 , alt_length )
  }else if( ref_length > alt_length ){
    ## DEL
    ref <- substr( tmp$Tumor_Seq_Allele1 , 2 , ref_length )
    pos <- pos + 1
    alt <- "-"
  }else{
    ref <- tmp$Tumor_Seq_Allele1
    alt <- tmp$Tumor_Seq_Allele2   
  }

  tmp <- paste( chr , pos , ref , alt , sep = ":" )
  result <- c( result , tmp )
}

dat_qc$POS <- result

##############################################################################

dat_maf <- dat_maf[,c("Chromosome" , "Start_Position" , "Tumor_Seq_Allele1" , "Tumor_Seq_Allele2" , "Hugo_Symbol" , "Variant_Classification")]
dat_maf$POS <- paste(dat_maf$Chromosome , dat_maf$Start_Position , dat_maf$Tumor_Seq_Allele1 , dat_maf$Tumor_Seq_Allele2 , sep = ":")
dat_maf <- dat_maf[,c("POS" , "Hugo_Symbol" , "Variant_Classification")]

dat_merge <- merge(dat_maf , dat_qc , by = "POS")
raw_mut <- nrow(dat_merge)

##############################################################################
## 1、突变的reads均不满足质控
if( length(which(is.na(dat_merge$Tumor_counts))) > 0 ){
  result_qc <- dat_merge[ -which(is.na(dat_merge$Tumor_counts)) ,  ]
}else{
  result_qc <- dat_merge
}

na_qc <- nrow(dat_merge) - nrow(result_qc)

##########################################
## 1、深度质控
dat_merge <- result_qc
result_qc <- dat_merge[dat_merge$Tumor_counts >= dat_stand$minimumdepth & dat_merge$Normal_counts >= dat_stand$minimumdepth,]
depth_qc <- nrow(dat_merge) - nrow(result_qc)

## 2、Tumor中突变reads的数量
dat_merge <- result_qc
result_qc <- dat_merge[dat_merge$Tumor_alt_counts >= dat_stand$minaltcount ,]
tumorAltReads_qc <- nrow(dat_merge) - nrow(result_qc)

## 3、Normal中突变reads数量的质控
dat_merge <- result_qc
result_qc <- dat_merge[dat_merge$Normal_alt_counts <= dat_stand$maxaltcount ,]
normalAltReads_qc <- nrow(dat_merge) - nrow(result_qc)


##########################################

##########################################
## 4、距离质控
dat_merge <- result_qc
result_qc <- dat_merge[dat_merge$DistanceToAlignmentEndMedian >= dat_stand$enddistance ,]
distance_qc <- nrow(dat_merge) - nrow(result_qc)

##########################################

##########################################
## 5、链偏倚质控
dat_merge <- result_qc
result_qc <- dat_merge[dat_merge$strandbiasprop <= dat_stand$strandbiasprop ,]
strand_qc <- nrow(dat_merge) - nrow(result_qc)
##########################################

##########################################
## 6、比对质量质控

## 1、0比对reads的质控
dat_merge <- result_qc
result_qc <- dat_merge[dat_merge$zeroproportion_Tumor <= dat_stand$zeroproportion & dat_merge$zeroproportion_Normal <= dat_stand$zeroproportion ,]
zeroproportion_qc <- nrow(dat_merge) - nrow(result_qc)

## 2、Tumor的reads比对质量
dat_merge <- result_qc
result_qc <- dat_merge[dat_merge$medianmapquality_altTumor >= dat_stand$minmapquality ,]
minmap_tumorAlt_qc <- nrow(dat_merge) - nrow(result_qc)

## 3、Tumor的alt的比对质量和Normal的ref的比对质量的差值
dat_merge <- result_qc
result_qc <- dat_merge[dat_merge$minmapqualitydifference_TumorAltNormalRef <= dat_stand$minmapqualitydifference  ,]
minmap_diff_qc <- nrow(dat_merge) - nrow(result_qc)

##########################################

##########################################
## 7、碱基质量质控
dat_merge <- result_qc
result_qc <- dat_merge[dat_merge$median_quality_Tumor_alt >= dat_stand$minbasequality ,]
minbase_tumorAlt_qc <- nrow(dat_merge) - nrow(result_qc)

## 存在突变为纯和的位点
dat_merge <- result_qc
dat_merge$median_quality_Tumor_ref <- ifelse( dat_merge$Tumor_alt_counts==dat_merge$Tumor_counts , dat_merge$median_quality_Tumor_alt , dat_merge$median_quality_Tumor_ref )

result_qc <- dat_merge[dat_merge$median_quality_Tumor_ref >= dat_stand$minbasequality ,]
minbase_tumorRef_qc <- nrow(dat_merge) - nrow(result_qc)

dat_merge <- result_qc
result_qc <- dat_merge[dat_merge$median_quality_Normal_ref >= dat_stand$minbasequality ,]
minbase_NormalRef_qc <- nrow(dat_merge) - nrow(result_qc)

##########################################
## 8、突变长度质控
dat_merge <- result_qc
result_qc <- dat_merge[dat_merge$mut_length <= dat_stand$mutlength ,]
mutlength_qc <- nrow(dat_merge) - nrow(result_qc)


##########################################
Variant_Type <- c("Missense_Mutation","Nonsense_Mutation","Frame_Shift_Ins","Frame_Shift_Del","In_Frame_Ins","In_Frame_Del","Splice_Site","Nonstop_Mutation")

## 9、质控SNV
# if candidate non-coding SNVs were in a tandem repeat region suggested by tandem repeat finder, we discarded the SNVs
dat_merge <- result_qc
result_qc <- subset( dat_merge , !( mut_length==0 & !(Variant_Classification %in% Variant_Type) & tandem_repeat_finder=="trf") )
SNV_NoncodingTRF_qc <- nrow(dat_merge) - nrow(result_qc)

# if candidate SNVs were in RepeatMasker repeat regions (http://www.repeatmasker.org) within 1Mb from the centoromeric or telemeric gaps, 
# we discarded the SNVs
dat_merge <- result_qc
result_qc <- subset( dat_merge , !( mut_length==0 & repeatmasker == "Simple_repeat" & telomere == "telomere") )
SNV_repeatmasker_telomere_qc <- nrow(dat_merge) - nrow(result_qc)

dat_merge <- result_qc
result_qc <- subset( dat_merge , !( mut_length==0 & repeatmasker == "Simple_repeat" & centromeres_1MB == "centromeres_1MB") )
SNV_repeatmasker_centromeres_1MB_qc <- nrow(dat_merge) - nrow(result_qc)

##########################################
## 10、质控INDEL
# if candidate non-coding indels were in repeat regions suggested by tandem repeat finder or RepeatMasker, we discarded the indels.
dat_merge <- result_qc
result_qc <- subset( dat_merge , !( mut_length!=0 & !(Variant_Classification %in% Variant_Type) & (tandem_repeat_finder=="trf" | repeatmasker == "Simple_repeat" )) )
INDEL_non_coding_qc <- nrow(dat_merge) - nrow(result_qc)

##########################################
## 11、性染色体质控
dat_merge <- result_qc
result_qc <- subset( dat_merge , !(Chromosome %in% c("chrX" , "chrY" , "X" , "Y") ))
Sex_qc <- nrow(dat_merge) - nrow(result_qc)

##########################################
dat_merge <- result_qc

dat <- data.frame( Raw_muts = raw_mut , 
  No_Alt_Reads = na_qc , 
  Depth = depth_qc , TumorAltDepth = tumorAltReads_qc , NormalAltDepth = normalAltReads_qc  ,
  DistanceToAlignmentEndMedian = distance_qc ,
  strandbiasprop = strand_qc ,
  zeroproportion_qc = zeroproportion_qc , minmapquality_TumorAlt = minmap_tumorAlt_qc , minmapqualitydifference_TumorAltNormalRef = minmap_diff_qc ,
  median_quality_Tumor_alt = minbase_tumorAlt_qc , median_quality_Tumor_ref = minbase_tumorRef_qc , median_quality_Normal_ref = minbase_NormalRef_qc ,
  mutlength_qc = mutlength_qc ,
  SNV_NoncodingTRF = SNV_NoncodingTRF_qc , SNV_repeatmasker_telomere = SNV_repeatmasker_telomere_qc , SNV_repeatmasker_centromeres_1MB = SNV_repeatmasker_centromeres_1MB_qc ,
  INDEL_Noncoding_TRF_Reatmasker = INDEL_non_coding_qc , Sex_qc = Sex_qc ,
  Muts_keep = nrow(dat_merge)
)

##########################################
write.csv( dat , out_record_file , row.names = F , quote = F )
write.csv( dat_merge , out_file , row.names = F , quote = F )